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Abstract. Dynamical topological solitons are studied in classical two-dimensional Heisenberg easy-axis 
ferromagnets. The properties of such solitons are treated both analytically in the continuum limit and 
numerically by spin dynamics simulations of the discrete system. Excitation of internal mode causes orbital 
motion. This is confirmed by simulations. 

PACS. 75.10.Hk Classical spin models - 75.30.Ds Spin waves - 05.45.-a Nonlinear dynamics and nonlinear 
dynamical systems 

1 Introduction act analytical solutions for 2D topological solitons in the 

isotropic magnet in the continuum limit, and proved that 



The analysis of two dimensional (2D) magnetic solitons 

S3: n 

, continues for more than 25 years, for reviews see Refs. Ill 

1 



Such solitons are well-known to play an important 
role in the physics of 2D magnetic systems. In easy-plane 
magnets with continuously degenerated ground state there 
appear magnetic vortices, which are responsible for the 
Berezinskii-Kosterlitz-Thouless phase transition Be- 
lavin and Polyakov were the first who constructed ex- 



such solitons are responsible for the destruction of the 

n 

long-range ordering for finite temperature jfj . In the ani- 
sotropic magnets such static solitons are unstable against 



collapse 



| . However in easy-axis magnets there appear 



various types of dynamical localized topological solitons 
due to the presence of additional integrals of motion. We 
will consider precessional solitons Q, which exist in uni- 
axial magnets due to the conservation of the z-projection 
Send offprint requests to: Denis D. Sheka 
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for a number of models used in field theory and condensed works of Papanicolaou et al |lJ . lfil lf| the dynamics of 



matter physics, see Ref. [lCj. The topological small radius 2D solitons was described using the algebra for some non- 
solitons become interesting now due to possible applica- canonic momentum. In particular, as it was mentioned in 
tions in high-energy physics and the quantum Hall Ref. ^| , a single topological soliton can not move without 

same time a free rotational motion 



m, , , c , , , . c , , . , ,., of the 2D topological soliton was predicted in Ref. 
I he problem of the dynamics of topological solitons ° 



17LIt 



, . ,. , , , , r n-rs r results in a finite mass for the small-radius soliton 17 1 , 

and vortices is a complicated task for 2D ferromagnets, *~ J ' 

i T , n . , mi while the mass of the localized soliton diverges as the log- 

where Lorentz and Galilean mvanance are absent. Ihe r— i ° ° 

r j. ,., arithm of the system size according to 181. 

presence ot a gyroforce acting on a moving soliton is the 1 — ' 

i j-u- i » i • , , , , , £ The present work is devoted to the analysis of the dy- 

only thing which is well established, but the free gyro- J 

j i , , , , .„ r r.^ namical properties of topological solitons in easy-axis fer- 

scopic dynamics nas not been reported till now tor any 2D 1 1 1 ° 

,., t-, ,-, , ill-, romagnets, both in the discrete model with weak anisotropy 

solitons. lor the easy-plane magnets the weak localization ° ' iJ 

r , , , . , , , and in the continuum model. We should stress here that 
ol tne vortex is related to tne gapless magnon spectrum; 

t , i j • i ,i i, ,i topological solitons were studied only in the frameworks of 

nence tne vortex dynamics is governed mostly by tne inter- 1 ° J 

.., , , , , , . , . , , . continuum field approaches in all above mentioned papers, 

action with the system border, and the mei'tial properties 

, 4 , .. j. Spin dynamics simulations for the motionless 2D topo- 

do not appear. As a result, computer simulations of mag- J 

, . . , , , c , ... , logical solitons were performed in 191. In this paper we 
netic vortex dynamics m a large but hmte lattice show ° !— " 

... , n -iii l-ii perform spin dynamics simulations for a wide range of 

a motion which can be described by complicated non J ° 

n 

»t • j . , , . .,, 1 Ji soliton shapes: from large radius solitons to small ones. 

INewtoman dynamical equations with nonlocal terms |lij. 

t-, ,i ..v ,i i j- The main issue is to move the soliton. Using the structure 

lor the isotropic magnet with the gapless magnon dis- i — i ° 

, iii j_i i_ r of internal modes 201, we have found such perturbations 

persion law numerical analysis shows the absence of the 

, n . i .. r , , ,., -, i of initial centrosymmetric soliton shape, which results in 
localized motion of the soliton |1J|. J r ' 

perfect orbital motion of the soliton. 

A very attractive candidate to discuss the general prob- 
lems of the magnetic soliton motion is the easy-axis fer- 

2 Discrete Model and Continuum Limit for 

romagnet. In this case the soliton shape is exponentially 

localized, it seems to be possible to separate the soliton 2D Ferromagnets 

motion from the magnons due to their finite activation TTT . , ,. . , , , „ , . , . , r 

We consider the simplest model of the classical 2D terro- 

energy. The general features of the 2D soliton dynamics, , , ., , , ,, <■ „ . TT -i. 

magnet described by the following Hamiltoman 

which should have particle-like properties with a finite T 

j#r = --J2( s n-Sn +a + ss* n si l+a ). (i) 

soliton mass, are not clear at present. For example, in (n,a) 
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Here S n = (S%, S^, S^) is a classical spin vector with fixed version of the Hamiltonian is 



length S (in units of the Plank constant H) on the site n of jq 2 r „ o sin 2 # 

*% <t>] = \ Jd 2 x (W) 2 + sin 2 9 (V0) 2 + ^ . (4) 

a two-dimensional square lattice, a is a vector to a near- L . 

i • 11 rp U i i • i j ,i ■ , • tt • u In terms of the fields 9 and d>, the Landau-Lifshitz equa- 
est neighbor. I he model includes the isotropic Heisenberg r ' H 

exchange interaction, J > is the exchange integral, and tions J5J) read 

«2 n 2 

(5) 



the spatially homogeneous uniaxial exchange anisotropy, 



8 is the anisotropy constant. The summation runs over 
nearest-neighbor pairs (n, n + a). The case (5 = corre- 
sponds to the isotropic model. To describe the anisotropy 
effects we will consider the case when 8 > 0, then the 
z-axis suppose the easiest magnetization. ui(k) — uq(1 + fc 2 Zg), (6) 



ftS 59 ' fiS , 

In the longwavelength limit the magnon excitations of the 
form 9 = const < 1, <j> = k ■ r — u>t, have the following 
dispersion law, 



The spin dynamics is described by the discrete version 
of the Landau-Lifshitz equations 
dS n 1 



dt h 



8. x °* 



BSr, 



(2) 



The model of the pure uniaxial ferromagnet has well- 
known linear excitations (magnons) above the ground state 
S z n = 1 of the form 1 - S z n = const < 1, S* + iS% oc 
exp(ik x a + ikya — iu)t), which have the finite gap disper- 
sion law 



which follows from J2J in the lowest approximation in 
ka < 1. 

3 The Structure of Precessional Soliton 



US 

u{k) = ujq H 7— 

n 



" n ) +Sm { — 



For the pure uniaxial ferromagnet the Hamiltonian |JI|I 
does not depend explicitly on the variable 4> due to the 
spin-space isotropy (in contrast to the lattice, which is al- 
ways anisotropic but in coordinate space). This condition 
(3) corresponds to the additional integral of motion 



Here ujq — 4JS5/H is the homogeneous ferromagnetic res- N — (S — S^) ■ (7) 

n 

onancc frequency, k is the wave vector. 



In the case of weak anisotropy, S <C 1, the character- 
istic size Iq — a/v4S of the excitations is larger than the 
lattice constant a, so that in the lowest approximation in 
the small parameter a/lo and with weak gradients of mag- 
netization one can use the continuum approximation for 
the Hamiltonian (TJJ by introducing the normalized spin 

S f 

s = S/S = (sin 9 cos fa sin 9 sin 0; cos 9) . The continuum N = — j d 2 x (1 - cos9) . (8) 



When N ^> 1 and the WKB approach is valid, one can 
consider N £ N as the number of magnons, bound in the 
soliton, see Ref. yj| . The conservation law (J7J) can provide a 
conditional minimum of the Hamiltonian, which stabilize 
the possible soliton solution, see below. The continuous 
version of JTJ) reads 
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The simplest nonlinear excitation of the model J^J is by p 2 d9o/dp and integrating over all p, one can easily ob- 



a 2D soliton, which has a finite energy. The topological tain the identity 
properties of the soliton are determined by the mapping 
of the xy-plane to the £— sphere of the order parameter 
space. This mapping is described by the homotopic group which gives a possibility to rewrite the soliton energy Q 



sin 2 9 (p)pdp ■■ 



2ujp 
Jo 



[l-cos9 Q (p)}pdp, (13) 



Tt2{S 2 ) = Z, which is characterized by the topological in- as follows 



variant (Pontryagin index) 



S ■ (V a S X V/3S) 



$ = (d C xc + flUJpN, 

JS 2 



(9) 



d 2 a; 



(14) 



The Pontryagin index takes integer values, q € Z, being 
an integral of motion. 

Let us consider the so-called centro symmetric topolog- 
ical precessional soliton, which has the following structure: 

= 0o (p) , = </?o + qx - Vpt, (10) 

where p = t/Iq is the dimensionless radius and uj p G 
(0, ujo) is the frequency of the internal precession. We will 
discuss only the case q = 1, when the soliton has a lower 



(V6>) + sin 2 9 (V4>Y 

Note that the linear dependence of the soliton energy § 
on N agrees with the general relation fku p = b~<S /SN. 

The shape of the soliton essentially depends on the 
number N of bound magnons. In the case of solitons with 
large radius R (equivalent to N 3> N ), the approximate 
"domain wall" solution works well, see Ref. jlj. This so- 
lution has the shape of a curved ID domain wall at the 
distance R 

r - R 

cos9 (r) = tanh — - — . (15) 



encr 



m, e r ,, r ,. n / \ i- n n r n Using this simple structure one can obtain the number of 
gy. I he form of the function #o(»J satisfies the follow- ° 1 



ing differential problem: 
d 2 6> 1 d9 



dp 2 p dp 



9 (0) - 7T, 



sm f o cos ( 



9 (oo) = 0. 



1 + - 

P 



-^-sin9 = 0, 
(11a) 



bound magnons, which is proportional to the area of the 
soliton, N w Nq(R/Iq) 2 , and precession frequency 



uq R 



(16) 



This equation was solved numerically in Refs. [21 
For the case of a centrosymmetric soliton the number of 
bound magnons 



(lib) In the case of small radius solitons (N <C Nq), the follow- 

n 

ing asymptotically exact solution works well 23] 



tan 



0o (r) _ R 



= - K x I - 



'•() 



= , (17) 



N = N Q pdp[l-cos9 (p)}, 



(12) 



JOj LU p /wo 

where K\{») is the McDonald function. It provides cor- 
rect behavior for r < R <Si Iq, where it converts to the 



where -/Vq = 2ixSl\ja is the characteristic number of Belavin-Polyakov solution tan6> /2 = R/r and provides 



3 



bound magnons for 2D magnets |21j . Multiplying Eq. (jllaf) a correct exponential decay for r > fl. In this case the 
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frequency of the soliton precession ui p — > loq when N — > 0, We will use this simple expression as initial condition 
but the dependence ui p (N) has a singularity at the origin: for our numerical simulations in Sec. [S] 
dajp/ dN — > co as N — > 

1 / 181 

oj ~ In (8iV /e7 2 7V) ' v ' 4 The Soliton Dynamics 



where 7 w 1.78 is the Euler constant, see Refs. 



3 3 



In the intermediate case of arbitrary R, it is possible To describe the dynamics of the soliton as a whole, it is 
to use an approximate trial function of the form, proposed necessary first of all to introduce an effective soliton coor- 



f. 



X(t) = [ d 2 x r(l-cos0). (22) 

Na z J 

Using this quantity we will look for the way of possible 
soliton deformation, which initialize its motion. In order 
to realize this idea, let us derive the soliton speed (see the 
Appendix A): 



Ref. 1211, dinate. Let us define the soliton position X(t) — X(t) + 

i\ / \ n / nX iY(t) as the center of mass of the z-component of the 

Oo(r) R I r — R\ ,„x 
tan = - exp — ) . (19) 

^ r V r o / magnetization field: 

n 

Here R is the fitting parameter, which was found in Ref. ;2Q] 
by fitting the trial function (|19|l to the numerical solution 
of the differential problem 111(1 . The value of this fitting 
parameter is closed to the soliton radius, which satisfies 
the condition cos6*o(i?) = 0. 

The trial function 1(19(1 gives a possibility to describe 
approximately the soliton shape for a given radius R. How- 

dX JS 2 ( 

ever, it contains one extra parameter, ui p , due to the de- — ^- = / d 2 irsin 2 #V</>. (23) 
pendence ro = ro(uj p ). One can calculate approximately 

It is convenient to classify all possible perturbations of the 

the w p (i?)-dependence as follows 

soliton shape using a complete set of functions. We choose 
w p(^) ~ R + l ' the solution of the linearized problem, which provides a set 

which provides the correct asymptote ((16(1 for R^$> Iq and °^ P ar ^ a ^ waves - 

gives the limiting value lo p = for R < l Q . In the same Let us remind that the soliton in an easy-axial ferro- 

approach the typical size of the exponential tail of the magnet has a number of local magnon modes. The exis- 

soliton r m lovf+p; thus for the soliton shape we tence 01 local modes is Possible because of the gap in the 



have finally magnon spectrum as predicted in Ref. 



such modes 



9 (r) R ( r — R I — R — \ correspond to different types of soliton shape oscillations, 

tan w — exp ' 1 



\ ^° V-R + ^oy ^21) To describe the local modes one has to linearize the Landau- 
fo + X — ^pt- Lifshitz equations JSJ) on the soliton background as it was 
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n 

done in Ref. [2fJ. We use the partial-wave expansion pendix A for the details) 

6 = 6 (p) + ^ A ™ ("m + V m ) COS# m , 
m 



(24) 



<Po + X ~ u p t + ^2 



sin t>o 



X(t) — i?orbe %UJ Roib — 

C-i =4 



\A-x\ \C^\N u l 



2NLU-! 



(26) 



where $ m = m X -u m t and co m is the magnon frequency in Thus ' onl y the Perturbation with the symmetry of the 
the rotating frame. The radial functions u m (p) and v m (p) mode ™ = -1 can lead to a soliton motion as a whole, 
satisfy the following eigenvalue problem for two coupled As we have found < the best wa ^ to excite such a mode is 



Schrodinger-like equations: 



d^__ l_d 
dp 2 p dp 



1 d 



dp z p dp ojq 



where the radial "potentials" are 



u m = W(p)v m , 
v m = W(p)u m , 



Tr . . (cos 6*n ± m) 2 
V± (p) = — + cos I cos V 



CO 



W(p) 



1 



to use the exact shape of this mode, calculated in linear 
approximation, with finite amplitude of deformation A-i. 
(25) We will check this prediction in Sec. El 

At the end of the section let us discuss the connection 
between the dynamics of the soliton center X(t), which 
results in the orbital motion l|26|l , and the dynamics of the 
specific soliton position, introduced by Papanicolaou and 



sin 2 O 1 



1 



2 



3 



?o) 2 - 



The local modes exist in a range of frequencies inside the 



Tomaras |15| as a some integral of the topological density 
Q ©: 

*=W- (27) 



gap, w m c G (0, ujQ — ujp). The number of local modes essen- This quantity can be interpreted as a "guiding center" of 
tially depends on the soliton radius: when the soliton ra- the soliton orbit. In the linear on A m approximation one 



dius decreases, the local modes leave the gap range, trans- can calculate (see the Appendix A) that 
forming to the quasi-local modes with singularities in the 

R(t) = —Ailoe x = const, 
scattering picture. For the soliton with R < R c = 1.52Zo 

there exists only one local mode, namely the mode with which corresponds to the simple soliton shift due to the 



(28) 



m = —1, and it is this mode which corresponds to the translational mode with m = +1. The physical picture 



soliton motion. 



is similar to the electron motion in the magnetic field: 



Let us calculate the effective soliton trajectory X(t) the electron moves along the circular Larmor orbit and 

using the partial wave ansatz (|24|l . In the linear approxi- dX / dt is not conserved; the generalized momentum P 

mation in A m all modes with \m\ 1 give no contribution also changes; however their combination, which determine 

to the integral (|23|l due to the angular symmetry, and the the center of the "guiding center" of the orbit R saves its 



effective soliton coordinat results as follows (see the Ap- position. 
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T3 

c3 



o 



o 

00 




We start the simulations using an initial soliton-like 
distribution 



20 40 60 80 

Number of bound magnons, N/N Q 



100 



Fig. 1. Radius of the soliton as a function of the number of 
bound magnons. The symbols correspond to the simulation 
data; the solid line is the result of the numerical integration of 
the continuum Eqs. 1111 . 11211 : the dashed line corresponds to 
Eq. ED 

5 Numerical Simulations for the Circular 
Symmetric Topological Precessional Soliton 



To validate predictions of the continuum theory for the 
soliton properties, we integrate numerically the discrete 
Landau-Lifshitz equations J2Jl over square lattices of size 
L x L using a 4th-order Runge-Kutta scheme with time 
step O.Of and periodic boundary conditions. In all cases 
the soliton is started near the center of the domain. We 
have fixed the exchange constant J = h = 1 as well as 
the spin length S = 1 . We have considered the anisotropy 
parameter in the range 5 € (0.0005; 0.1), corresponding 
to lo/a £ (22.4, 1.58) so that we are close to the contin- 
uum limit. We consider system sizes in the range L/a G 
(50,800). 



9 (r), c/) = ip +x, 



(29) 



with the trial function 121fl for the ^-field. Evidently, the 
soliton solutions of the Landau-Lifshitz equations for a 
lattice differ from the circular symmetric continuous solu- 
tions, and also from the simple trial function. To find such 
a "pure" soliton solution, i.e. to adapt the trial solution to 
the lattice, one should provide enough time for the decay 
of the initial error in the trial functions. In fact, using <|29[l 
as initial conditions for the lattice we excite also magnons, 
which should be taken out the system. To avoid the prob- 
lem of magnons we have damped them the initial stage of 
simulations by applying damping. This kills all spreading 
spin waves coming from the imperfect initial condition. In 
this way instead of Eqs. (J2J), we have integrated numeri- 
cally Landau-Lifshitz equations with Gilbert damping 



\S x 


e 




\ S x 




[ n dS n _ 


+ S 




n dS n _ 





see Ref. 25] for details. The lowest frequency of the con- 
tinuous magnon spectrum is loq, thus the damping time 
id ~ l/( etJ o)j see R°f- 19] for details. During the damp- 
ing time (t < td) , the magnons are damped in the system, 
but the soliton is also damped, and the soliton energy S 
decays as well as the number of bound magnons N. In or- 
der to save the soliton structure, we should switch off the 
damping before we damp out the soliton, i.e. t < l/eu> p . 
In all simulations we use the same value of e — 0.02, then 
the damping time td w 12/5, and the damping is turned 
off adiabatically after a time greater than td- 
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(a) 



(b) 



o 

C/5 



100 200 300 400 

Number of bound magnons, N/N Q 

Fig. 2. Soliton precession frequency as a function of the num- 
ber of bound magnons. Symbols correspond to the simulation 
data; the solid line is the result of continuum model integra- 
tion; the dashed line corresponds to Eq. 132II . 

Let us discuss the choice of the other parameters. In 
all simulations we want to be not far from the contin- 
uum limit in order to validate the continuum approach. 
It means that the magnetic length l should be greater 
than the lattice constant a. This regulates the choice of 
the anisotropy constant S = a 2 /4Zq. Besides Iq the soliton 
shape is characterized by two extra scales: R, which is the 
soliton radius, and ro, which characterizes the scale of the 
exponential decay of the excitation far from the soliton 
center. 

We start with the large radius solitons. In this case 
R S> ro w Iq, and we can limit ourselves by choosing 
5 = 0.1 (this corresponds to Iq « 1.6a). The system size L 
should be much greater than the largest parameter of the 
soliton, which is its radius. We consider solitons up to the 
radius R — 201q ~ 31.6a. Thus we consider lattices with 
L = 200a , which satisfy all above mentioned conditions. 



. «• *• *■ 

tfrf » 

4 ' 



% 



' 2> 



\ < 




20 22 24 26 28 30 32 34 



20 22 24 26 28 30 32 34 



Fig. 3. In-plane spin distribution for the soliton with R = 5lo, 
which is situated at Xo = 26.3 + i25.0. Fig. (a) corresponds to 
the circular symmetric soliton and (b) to an elliptically de- 
formed one. The lines describe the contour plot of S z — 
(6 = vr/2): the dashed line is for the circular symmetric soliton 
and the solid line is for the deformed soliton. 

In the case of small radius solitons we have the follow- 
ing relation between the parameters of the system: 



a < R < Z < r < L. 



(30) 



For the smallest soliton we choose <5 = 0.0005, which corre- 
sponds to Iq w 22.4a; this gives the possibility to consider 
solitons of small radii down to R = 0.225/o ~ 5a. How- 
ever, such a small anisotropy drastically changes the soli- 
ton shape far from the center, which has the scale ro, see 
Eq. ((TT}l . For example, for the soliton with R = 0.225Z , 
the precession frequency cu p ph 0.84wo (see Ref. [2fj), which 
results in ro = lo/y/1 — uj p /ujq ~ 56a. Thus to consider 
small radius solitons we must increase the system size. 
In our spin dynamics simulations we choose L = 800a 
for the smallest solitons. To perform simulations for such 
large systems, 800 x 800, we have used parallelize a com- 
putations, see the Appendix B. 
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Let us discuss results of our spin dynamics simula- 
tions. Starting from initial conditions (|29|l . and adapting 
the soliton shape to the lattice, we have obtained the class 
of one-parameter stable soliton solutions for a wide range 
of the parameter N, or equivalently, the soliton radius R. 
We have studied the R(N) dependence, which is presented 
in Fig.^ Almost in the full range of parameters, the sim- 
ple dependence 



N ^ N Q 



(31) 



is valid. Note that the dependence (|31[) was verified nu- 
merically in Ref. for the large radius solitons only, 
where R > lOZo- Here we want to check the continuum re- 
sults for arbitrary R. Using (|2"0)l and (|3~T]l , an approximate 
dependence 



(32) 



1 + ^N/No 

can be derived. To compute the precession frequency, we 
calculate the Fourier spectrum of the in-plane spin compo- 
nents. One can see from Fig.|21that this simple dependence 
works well in a wide range of parameters. 

6 Simulation of the Orbital Motion of the 
Soliton 

In the previous section we have performed spin dynamics 
simulations only for the circular symmetric precessional 
soliton, which does not move as a whole. As mentioned in 
Sec. rj] in order to move the soliton one should break its 
symmetry. We have done this by an initial deformation of 
the soliton shape, and integrated numerically the Landau- 
Lifshitz equations. Specifically, we have chosen an elliptic 



26 



>h 25 



24 




24 



25 
X 



26 



Fig. 4. Trajectory of the soliton with initial radius R — 4.87Zo 
and initial precession frequency u) p — 0.2ajo- The initial point 
Xo+iYo corresponds to the center of the soliton from the Fig. [3] 

kind of deformation, which corresponds to the shape of 
the internal partial mode with azimuthal quantum number 
m = — 1. We start the simulations with initial conditions 



9 (p) +A[u- 1 {p) +v-i(p)] cos X, 
A 



(33) 



4> = <A) + x - 



sin t>o 



[u-i(p) - i;_i(p)]sinx 



by the same numerical scheme as described in the pre- 
vious section. We calculated the functions U-i(p) and 
V-i(p) numerically solving the eigenvalue problem (|25|l 
by the two-parametric shooting scheme as described in 
Ref. 2^- 

The parameter A is the amplitude of the eigen- 
mode, which characterizes the magnitude of the soliton 
deformation. An initial distribution of spins, which corre- 
sponds to Eqs. I|33|) . is shown in Fig. [3] and can be seen to 
describe the elliptical kind of the soliton deformation. 



10 
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0.2 0.4 0.6 0.8 1 

Soliton precession frequency, G) p /G) 

Fig. 5. Frequency of the orbit motion as a function of the 
precession frequency of the soliton. Symbols correspond to the 
spin dynamics simulation data; lines corresponds to the fre- 
quency of the eigenmode with m = — 1 from the continuum 
theory by Sheka et al. 



uue \ 

Q 



During the simulations we have computed the time 
dependence of the position X(t) — (X(t),Y(t)) of the 
soliton center: 



X(t) 



(34) 



Tim & Sn(t) 

which is the discrete analogue of Eq. (|22[1 : r n — (x n ,y n ) 
are the lattice points. 

We have found numerically that after switching off the 
damping, the soliton reaches very fast a circular trajec- 
tory, see Fig. 01 This results in a nice circular motion with 
constant frequency. For small deformations the radius of 
the orbit is proportional to the initial deformation. One 
can say that the effect of a circular motion and the excita- 
tion of the mode with m = — 1 are identical for this case, 
as predicted by Sheka et al. [^(J. Such a relation is valid 
in some range of the soliton deformation for all soliton 




1 2 3 

Amplitude of deformation, A 

Fig. 6. Radius of the soliton orbit as a function of the defor- 
mation amplitude A. Parameters of the soliton: initial radius 
R — 4.87/o and initial precession frequency ui p — 0.2uo. 

radii, see Fig.[5J Then, for larger deformations, non-linear 
regime is clearly seen, see Fig. El The frequency of this 
orbit motion of the soliton approximately corresponds to 
the frequency of the local mode, f2 rb = w m= _i. 

The presence of such an exact circular motion, with 
only one frequency, independent of the orbit radius (even 
in non-linear regime), gives the possibility to conclude 
that this is the first observation in the numerical experi- 
ment of the pure gyroscopic motion, which is equivalent to 
the Larmor precession of a charged particle in a magnetic 
field. Therefore, the soliton motion can be described by an 
effective equation of motion for the position of the soliton 
X, which takes the form of usual Newtonian equation for 
a particle with the well-defined effective mass M under 
the influence of the gyroscopic force F g : 



F g =G 



e z x 



dX 



(35) 
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Here G = AirHS/a 2 is the gyroconstant, see Formally, mass increases faster that the domain wall width, which is 

Eq. H35|) has two solutions. One solution, X = const, cor- proportional to R. The mass takes a minimum value M c « 

responds to the translation mode with m = +1. In the 3.51Afo for the soliton with R c ss 0.547Zo The soliton with 

infinite system this is a zero-frequency local mode, which these parameters has the highest mobility. 

describes a simple shift of the soliton, ui m=+ i = 0. The 

second solution describes a circular motion with the fre- J Conclusion 

quency Q or h — G/M . Thus, we can calculate the effective 

In this paper we have studied the dynamics of topological 

mass of the soliton using the simulation data for the orbit 

solitons in classical 2D easy-axis ferromagnet. The analy- 

frequency. 

sis was made both analytically in the continuum approx- 

We have checked how the orbit frequency depends on 

imation and numerically using the spin dynamics simula- 

the soliton radius. For the large radius solitons i? rb ~ 

tions for a wide range of solitons: from large radius solitons 

2ojq(Io/R) 3 ' is in good agreement with the results for the 

to small ones. Our simulations were performed for small 

local modes. For the small radius solitons, the orbit fre- 

anisotropies, which corresponds to the continuum descrip- 

quency tends to the boundary of the spectrum, J7 rb ~ 

tion. We have checked and confirmed a number of results 

u>o — u)„. This dependence corresponds to our result for the 

Qfrom the continuum theory about the soliton structure, in 
. In the case 

particular, the connection between the number of bound 

N <C Nq one can use the approximate limiting expression 

magnons and the precession frequency of the spins inside 

<|18l) for uj p , which results in I2 rb ~ ^o/ m (8-/Vo/e7 2 iV). 

the soliton. 

The frequency tends to zero when R — > 0. 

The main issue is connected with the soliton dynam- 

One can calculate the effective mass of the soliton by 

ics. We have proposed a way how to move a soliton excit- 

M = G/(2 aih , 

ing one of its internal modes. To our knowledge, it is the 

N\ „ G nfi 2 



^ ~ ^IqF [ iVo/ ' ~ ujo ~ Ja^S 1 nrs t observation of inertial motion of 2D magnetic soli- 

where Mo is a characteristic value of the effective mass, tons. The effective soliton dynamics is similar to the Lar- 

found in Ref. 17]. The function F(») = ojq/ fi tb{*) de- mor dynamics of a charged particle in a magnetic field, 

pends on the soliton size having the asymptotic behavior By analysis of the effective soliton dynamics we extract 

information about the effective mass of the soliton. This 

{fi 1 I ) when mass essentially depends on the anisotropy, M oc 1/5, 

\ejxj (3?) » I 

1^,3/2 when j> 1 an( f on t ne soliton size, having the minimum for the soli- 

The soliton mass diverges in the limiting cases when R — > ton of the radius about 0.5/o- In the case of large radius 

and R — > oo. Note that for the large radius solitons the solitons the soliton mass increases with the increase of 



12 Sheka, Schuster, Ivanov, Mertens: Dynamics of topological solitons in two-dimensional ferromagnets 

the soliton radius. Note that it increases faster than the Appendix A Calculation of the soliton velocity 

number of the bound magnons, M ~ ^M (N/N ) 3 / 2 . 

In order to check the possible soliton motion we calculate 

Such dependence is in a good agreements with previous 

the speed of the effective soliton position X(t). According 

results [lj |; the soliton mass diverges in the limit case 

to the Eq. ll'L'l) the soliton speed 
R — > oo, which corresponds to the fact that the single S f If 5S 

~aT = 7v~~2 I ^ x Xi sin^fl = — d 2 x Xi — , 
domain wall can not move. When the soliton radius is at i\ ci J ni\ J o<p 

smaller, the effective mass increases logarithmically with where we used the Landau-Lifshitz Eq. ©. Calculating- 
its radius, and diverges in the limit case R = 0. Note that the functional derivative for the energy functional Q and 
the problem of inertial properties of a small radius soli- integrating by parts using the identity 
ton have caused a lot of discussions. According to linear Xi V • (sin 2 0V</>) = V • (a* sin 2 6V<f) - sin 2 93^, 
analysis 171 . the soliton mass tends to some limit value 

one can derive the soliton velocity in the form 
M* = 6y/ir(ir + 2)M w 55M when R — > 0. At the same , . „ 2 



(38) 



rH dX 2tt1 JS 2 f°° , / . 2 

time our previous analysis of cigcnmodes 20] shows that ar niy Jo v 

M -> oo at the limit case R -> 0. Spin dynamics Simula- which is equivalent to E q- 03 • Here the averaging means 



(F(.,x)) = (Wf^F(., X )dx 



tions confirms our results on internal modes: the soliton 



loose its mobility when becomes very small. Usin S the Partial-wave ansatz C2I, one can concretize 

the average value in Eq. 1381) : 



sm Wq 



sin <P m e lx 



We have predicted the fine circular motion of the soli- ^ s i n 2 OVpCpj = ^4 m <sin 2 9q 

m L 

ton by exciting its internal mode. We believe that such 
phenomenon can be observed experimentally, e.g. by ac 
pumping. Our investigations can be important also for the Aft er averaging with account of the expressions 
quantum Hall systems, where skyrmion-type solitons are (cos<? m e* x 



+ - [(% + v m ) sin20o + m(u m - v m ) sm9 ] ^cos^ m e 4X ^> 



2 

well-known to lead to the breakdown of the spin-polarized ( sm<P m e ix \ — \ m ^ 1 c imuj m t 



(39) 



quantum Hall effect [2Jj . one can calculate the soliton velocity in the form 

dX iluqIqNq ( i t -iu_ x t\ f An \ 

n = ^r\ AlCie + A -^ e J> ( 4 °) 

D. D. Sheka thanks the University of Bayreuth, where part wh ere the constants C m are determined by the static soli- 

of this work was performed, for kind hospitality and acknowl- ton structure, 

edges support from the Alexandr von Humboldt Foundation. Cm = _ 2 j pdp cos Q Q J (j + m ) (u m Ui - v m v^j 

C. Schuster thanks Stefan Karpitschka for his kind support ~. 

concerning aspects of parallelizing the source code. ^ ^( rnVl m Of 
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The mode with m = 1 is the zero-frequency local mode, that 
which describes a shift _of_the soliton position, its eigen 
spectrum has the form 



. U m V— m VmlA—m 



) 



LO {) 

V ■ (u m Vv- m + u m Vit_ TO - w- m Vf m - V- m Vu m 



ui=9 , v 1 =6 -\ , wi=0. 

P P 

for any m. The righthandside is in the form of the total 

A simple calculation shows that C\ = 0, therefore the ... .. . 

divergency, thus it gives no contribution to the integral 

soliton motion is connected only with the mode m = — 1, 

over the system: 

see Eq. (EH!) • x 

t , i i , M l « , „ / Pdp {u m v- m - v m u- m ) = 0. 

Let us calculate the dynamics of the guiding center J 

position of the soliton (gZJ. Using the partial-wave expan- Therefore the constant C*_ 1 =0, and finally the "guiding 

sion one can rewrite the topological density Q © as center" of the soliton can be rewriting in the form 

follows: 

Q = — smd 9 ° 9 '° - J- V A m cos $ m \(u' + v') sin O Appendix B Parallelized Computations 

l oP l oP"zf L 



+ (u m + v m ) cos 6*06*0 + m(u m - v m )b 
Averaging the linear momentum (|27|) with account of l|39fl 



There exist two main possibilities of parallelization com- 
putations: (i) the usage of a vector-computing machine 
which has a shared memory or (ii) the usage of a cluster 

one can derive R(t) in the form 

system (usually a Linux cluster) with communication be- 
R(t) = IqA-lCI + IqAxC^ exp (-iuj^t) , (41) tween the different p rocesS ors (nodes). The latter process 

where the constant can be calculated as follows: is called message passing interface (MPI), this interface 

does not depend on the programming language. 



C* = — 

m 4 



\j P d PS (1 + m){u m ui - v m v^j 

I In MPI there is a master-process which is responsi- 



- (1 - to) (u m vi - v m u^jj. 
One can easily see that C{ = —1, and the contribution 



blc for the administration of the data, i.e. initializations, 

reading or saving of data, whereas the other nodes (slaves) 

„ ,. , , „ . ,, , .„, „ ;1 are doing the calculation (e.g. integration). In our case we 

of the mode m = +1 results m the soliton shift, tor the 

... . . . , divide the lattice into horizontal stripes. As in our system 

mode with in = — 1 one can rewrite tne constant C_ 1 111 

. r we take into account nearest-neighbor interaction we must 

the form: 

put a communication between the borders of the stripes. 

C_i = — I pdpl U-\Vi — V-\Ui ) . 

2Jq\ I To make calculations for the i — th stripe we need the lower 

This integral vanishes due to the symmetry of the eigen- border of the (i + 1)— th stripe and the upper boarder of 

value problem 1)25(1 . Namely, a simple calculation shows the (i — 1)— th stripe. This exchange is done after every 



14 Sheka, Schuster, Ivanov, Mertens: Dynamics of topological solitons in two-dimensional ferromagnets 

second integration step in the Runge Kutta algorithm. The parallelizing of the source code is much easier 
This advantages a good equilibrium between the latency and more effective if one can use shared memory multi- 
period (time in which the nodes are synchronized), data- processors machines. The loops are environed by sunstyle 
transferring time, and calculation time. As the boundary parallelization-directives. According to the dependency of 
condition of our lattice is periodic the upper boarder of the variables they have to be declared as private variables 
the top stripe is exchanged with the lower border of the and others as reduction variables. Private variables take 
bottom stripe. Concerning the left and right boarders of different values in different threads (such as auxiliary vari- 
one stripe there are no communication processes, because ables) so there is no shared memory for these variables. Re- 
horizontal cuts of the system advantage an internal ex- duction variables are used for sums which are distributed 
change of the lateral borders in one process. over different threads during one loop and are added in the 

end of the parallelized loop. The advantage of this method 

As the data is always transferred between the same 

is a very fast access to the shared memory of all threads. 

nodes we implement a persistant connection mode in or- 

This avoids a time consuming latency period caused by a 

der to make the overhead (additional time spent on the 

communication process. To summarize we want to stress 

connection establishment) smaller. Theses connections are 

that a shared memory machine is always much faster than 

built up at the beginning of the integration to be ready 

a cluster, if one considers a fixed number of processors to 

for a fast use. 

be in use. The only advantage a cluster has is the fact that 
The transfer takes place in an asynchronous, buffered compared to a shared memory machine it has much more 
mode. Thereby the data resulting of the calculation on a processors and the number of processors can be increased 
stripe is put into a buffer so that for this node there is no arbitrarily. Also aspects of the price motivate universities 
need to wait until the other nodes are ready to receive. In and institutes to buy a cluster instead of a shared memory 
order to do a new calculation in the next time steps the machine, 
stripe awaits information of its neighbored stripes. The 
advantage of this method is that the time exposure for References 
the synchronization of the nodes decreases. 
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